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Abstract. We present exact analytical calculations of scanning tunneling currents in locally 
disordered graphene using a multimode description of the microscope tip. Analytical 
expressions for the local density of states (LDOS) are given for energies beyond the Dirac 
cone approximation. We show that the LDOS at the A and B sublattices of graphene are out 
of phase by tt implying that the averaged LDOS, as one moves away from the impurity, shows 
no trace of the 2qF (with qp the Fermi momentum) Friedel modulation. This means that a 
STM experiment lacking atomic resolution at the sublattice level will not be able of detecting 
the presence of the Friedel oscillations [this seems to be the case in the experiments reported in 
Phys. Rev. Lett. 101, 206802 (2008)]. The momentum maps of the LDOS for different types 
of impurities are given. In the case of the vacancy, 2qF features are seen in these maps. In all 
momentum space maps, K and K -\- K' features are seen. The K -\- K' features are different 
from what is seen around zero momentum. An interpretation for these features is given. The 
calculations reported here are valid for chemical substitution impurities, such as boron and 
nitrogen atoms, as well as for vacancies. It is shown that the density of states close to the 
impurity is very sensitive to type of disorder: diagonal, non-diagonal, or vacancies. In the case 
of weakly coupled (to the carbon atoms) impurities, the local density of states presents strong 
resonances at finite energies, which leads to steps in the scanning tunneling currents and to 
suppression of the Fano factor. 



New J. Phys. 
1. Introduction 

Graphene [[H O consists of a monolayer of covalently bonded carbon atoms forming a two- 
dimensional honeycomb lattice O lU. Low-energy electronic excitations in graphene are 
well described as massless Dirac fermions with an additional pseudospin degree of freedom. 
Because of the Dirac spectrum, impurities can have a strong effect on the local electronic 
structure of graphene when the Fermi energy is near the Dirac pointBSl [6l CI El 111 [IHl [TT1| . 
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Impurities in graphene can be in the substrate, in the form of adatoms, or as imperfections in 
the lattice itself. In one hand, there has been very significant progress in decreasing the amount 
of disorder introduced in graphene, for example by fabrication of suspended samples O. 
On the other hand, impurity effects have been explored to modify and tailor the electronic, 
thermal and chemical properties of graphene. Examples of the later include experiments with 
graphane[fT3l, chemical substitution of some of graphene's carbon atoms by boro! n an! d 
nitrogen atoms [I T4l[T5]l . and doping of graphene with metals on top [fT6l[T7l . 

Since graphene is an atomically thin membrane, it can be easily accessed with Scanning 
Tunneling Microscopy (STM) measurements. Impurity effects can be studied with atomic 
resolution and the local spectrum can be obtained by STM spectroscopy. In addition, atomic 
manipulation can also be performed with STM. There has been several STM studies of 
graphene grown epitaxially on SiC lfT9l [20l [2T]| . mechanically exfoliated graphene on Si02 
lEIl [H El [25], and graphene flakes on graphite[26J. In fact, STM experiments have 
proved instrumental in mapping the topography of corrugated graphene and determining the 
existence of charge puddles [l27l . Additionally, STM experiments are also able to probe 
the chiral nature of the electrons in graphene when they scatter from impurities [[28]|. This 
experimental work showed that intravalley backscattering is virtually absent in graphene. In 
particular, the STM experiment showed the lack of the 2qF (q^ is the Fermi momentum) 
Friedel modulation on the local density of states of graphene. As we show explicitly below, 
this lack of modulation can be traced to the fact that the local density of states at the A and 
B sublattices are out of phase by tt, an aspect already noted in passing previously [|29l[30l . 
Therefore, the local density of states, when averaged over the unit cell, shows no trace of 
the 2qF oscillation. Additionally, as we show below, at distances d close to the impurity, 
d <C l/qr^ there is a strong departure from the spatial dependence JlIlOl of the local 
density of states. Moreover the form of the density of states close to the impurity is very 
sensitive to type of disorder: diagonal, non-diagonal, or vacancies. 

In this work, we present calculations of STM currents in locally disordered graphene. 
We focus on the case of chemical substitution (by boron or nitrogen atoms, for example) 
and use a multimode description for the STM tip. We obtain exact analytical expressions 
for the local density of states, and also present results for energies beyond the Dirac cone 
approximation. We model the substitutional impurity by both an on-site impurity potential and 
local hopping disorder. We find that inclusion of the hopping disorder term leads to additional 
higher harmonics oscillations in real space for the density of states for the sub-lattice that 
does not contain the impurity. The main oscillations in the two sub-lattices are out-of-phase 
away from the impurity. For the regime in which the electronic hopping between the impurity 
and the nearest neighbor carbon atoms is decreased in relation to the hopping between carbon 
atoms in the clean system, the local density of states presents strong resonances. A vacancy 
is a extreme case of this regime. These resonances lead to the appearance of steps in the STM 
current, a signature that should be observable experimentally. These resonances also leads to 
open channels for tunneling between the STM tip and graphene, and lead to a decrease in the 
Fano factor. 

Another issue addressed in this paper relates to the effect of the tip on the measured STM 
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currents. In the usual analysis, the electrons in the tip are represented by jellium model with 
constant density of states. In this type of model neither the real part of the self-energy due 
to the tip-system coupling nor the variation of the density of states with energy is included 
(wide band limit). The tip, however, is not an infinite metal. In fact it has a structure where 
the number of atoms in the atomic planes reduces as we approach the tip. In a previous 
publication [3 2 J we have modeled the tip as a one dimensional model (in that work we have 
also considered the simplification of zero on-site energy at the impurity), which corresponds 
essentially to the case of a constant density of state too a good approximation. In that case we 
found the STM current to be symmetric around zero energy. When we generalize to the case 
of a multimode tip this symmetry is lost, as we show in this work. Comparing the results of 
Ref. [|32II with those given here it is possible to disentangle the effects due to graphene and to 
the impurities from those due to the tip. This work shows that some care has to be taken when 
interpreting the STM currents directly. 

The present manuscript is organized in the following way: In Section 2, the Green's 
function formalism is presented, with analytical results for the Green's function for graphene 
with a substitutional impurity, and for the STM tip modeled by a multimode system. Local 
density of states results are presented in Section 3, and results for the STM current are 
presented in Section 4. Section 5 contains final discussions and conclusions. 

2. Graphene and STM tip Green's functions 

2.7. Graphene 

The honeycomb lattice has two carbon atoms per unit cell, one from sublattice A and one 
from sublattice B, as depicted in Fig.[Il The unit cell vectors are ai and a2, with magnitudes 

= \a2\ = a, where a = ^/3ao 2.461 A, and ao is the carbon-carbon distance. Any 
lattice vector r can be represented in this basis as r = nai + ma2, with n, m integers. In 
Cartesian coordinates, ai = ao(3, \/3, 0)/2 and a2 = ao(3, —a/3, 0)/2, and the reciprocal 
lattice vectors are given by: bi = 2^(1, ^3, 0)/(3ao) and 62 = 27r(l, -V^, 0)/(3ao). 
The vectors connecting any A atom to its nearest neighbors are: 61 = (ai — 2a2)/3, 
62 = (a2 - 2ai)/3, and ^3 = (ai + a2)/3. 

We consider here the case where a substituting atom replaces a carbon atom in the A 
sublattice, say. When this happens two effects take place: (i) the on-site energy at the 
impurity site is different from that at the carbon atoms; (ii) the hopping from and to the 
impurity atom, t^, changes relatively to that of pristine graphene. In the latter case, we model 
the change in the hopping by introducing an additional non-diagonal term to the Hamiltonian, 
such that ti = —t + to (see below). Using these definitions the Hamiltonian can be written as: 
H = Ho + Vt + Vi, where 

Ho= - t^[b\r)a{r) + b\r - a2)a{r) + b\r - ai)a{r) + h.c], (1) 

r 

is the kinetic energy operator and a (6^ b) are fermion creation and annihilation operators 
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Figure 1. (color on line) Honeycomb lattice of graphene, with an substituting impurity at the 
A sublattice (the square). The unit cell vector ai and a2 as well as the next nearest neighbors 
vectors Si (i =1,2,3) are also represented. 

in the A (B) sites. The spin index is omitted for simphcity. We consider an isolated impurity 
located at r = (0, 0, 0), on sublattice A, so that its contribution to the Hamiltonian has two 
terms: 

Vt = to[b\0)a{0) + b\-a2)a{0) + b\-ai)a{0) + h.c] (2) 

and 

V, = e,a\0)a{0) , (3) 

for hopping and potential disorder, respectively. In the case of zero chemical potential, when 
the Fermi level crosses the Dirac point, the system is most susceptible to the presence of 
impurities. In the particular case to = hopping to the impurity site is completely suppressed, 
and the scattering term Vt represents a vacancy. It is well known [34J that the formation 
of a vacancy will lead to some local distortion of the carbon-carbon bonds. This effect is 
not incorporated in our Hamiltonian, which in that case would not be exactly solvable. The 
substitution of carbon atoms by boron or nitrogen has the main consequence of changing the 
local hopping and the onsite energy, both effects included in our description. The particular 
choice for boron or nitrogen is due to size restrictions imposed by the unit cell of graphene. We 
note here that replacement of carbon atoms by boron and nitrogen was already experimentally 
achieved L15.1 . 

We first calculate the single particle Green's functions for the system comprised of 
graphene and a single impurity, described by the Hamiltonian H above. The single particle 
Green's functions carry sub-lattice indices, and are defined as: 

Gaaik,q,T)= -(r[afc(r)at(0)]> , (4) 
Ghb{k,q,T) = -(r[6fc(r)6t(o)]> , (5) 
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Gab{k,q,T) = -(r[afc(T)6t(0)]> , 
GUk,q,T) = -(T [Mr) 4(0)]) . 
The equations of motion for the Green's functions are given by: 



Q 



iuJnGah{^n,k,p) = 



>^k,qGha{^n^q,P) + -^G aa{^n, P) 



E 



Xk,qGtb{(^n^ <7, P) + -^Gab{(^n^ 9, P) 



where 



Afc,p — —t(pp{Sk^p — to /Net) , 



(6) 
(7) 

(8) 
(9) 

(10) 
(11) 

(12) 
(13) 



and Nc is the total number of unit cells in the lattice, and Un are fermionic Matsubara 
frequencies. Note that Xpq ^ Xqp. This is a consequence of the impurity hopping term 
Vt, which breaks sub-lattice symmetry. The impurity potential term Vi also breaks sub-lattice 
symmetry, and therefore appears in an asymmetric way in the equations above. The set of 
equations of motions can be solved exactly. The presence of the scattering term Vt leads to the 
appearance of the phases 0fc and a more complex form for the T-matrix than usual. The exact 
solution for the Green's functions can be written, after a lengthy calculation, in the form LlOl : 

Gaaik,p) = Gl + g + h [Gl + Gl] + Gl T G^ , (14) 



G,,{k,p) =S,,pGl + ^-^GlTGl^-^ 



(15) 



where all the terms (G, G^, g, h and T) also depend on ujn (omitted here for brevity). The 
terms g, h, and T correspond to sums over infinite series of Feynman diagrams for impurity 
scattering, and are given by: 

g{ujn) =tlG\uJn)/[N,D{uJn)l (16) 
h[ujn) = to{t - to)/[N,D{u;n)l (17) 

and 

iUnto{2t - to) - Cit^ 



where the denominator D{ujn) is defined as 

D{ujn) = {t- to)' + [iuJnto{2t - to) - e,t^] 

and 



(18) 

(19) 
(20) 
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with the diagonal component of the Green's function for the clean system given by (h = 1) 

Gk = G^i^n, k) = 721X12 ' ^^^^ 

which is translationally invariant. The expressions for the Green's functions, Eqs. ([T4h-([T9l), 
are exact analytic solutions for graphene with one isolated substitutional impurity, including 
contributions from both the on-site energy and the hopping parameter to- Inclusion 
of the off-diagonal disorder to leads to additional terms, and additional cj-dependence of 
the graphene Green's function. Since single particle properties, such as local electronic 
spectra, can be obtained directly from the Green's functions, this cj-dependence has direct 
experimental consequences, such as for STM spectroscopy measurements. These results for 
the graphene Green's functions have been obtained in Ref. ifTOl . where local density of states 
maps, electronic spectra and Friedel oscillations have been calculated for the cases of boron 
and nitrogen substitution. We include here a brief derivation of the analytical expressions 
for the Green's functions, Eqs. ([T4l)-(l2T1). for completeness. Upon closer inspection, the 
physical meaning of the extra terms in Gaa become clear. The significance of the term g{(jJn) 
in ([T4l) which only appears in Gaa^ is more easily interpreted if we do a double Fourier 
transform to real space. This term corresponds to the return amplitude to the impurity site 
for an electron starting at the impurity site. The factor 1/ D{ijJn) contains a sum over an 
infinite series of intermediate scattering events, but the overall process is bounded and the 
tl factor denotes hopping from the impurity to the nearest neighbor 5-sites and back to the 
impurity site. Likewise, an interpretation can be given to the other term which only appears 
in Gaa^ namely, h{uJn)G^j^{(jJn)' A double Fourier transform shows that this term contributes to 
Gaai'f^i 0) and describes the amplitude of propagation between the impurity site and another 
A site, again with an infinite series of intermediate scatterings. Similarly, the h{uJn)Gp{uJn) 
term contributes to Gaa(0, r). No such terms can, of course, appear in G55 when the inpurity 
is at a ^4 site. And no such term can be present when there is only the impurity potential term 
Ci. The G^{uJn)T{uJn)Gp{uJn) term, which appears in both Gaa and is the usual term also 
present in simple on-site impurity potential problems, but in this case the T-matrix contains 
contributions from both to and e^. 



2.2. STM Tip 

Let us consider a model for the STM tip represented by a multimode system. The bulk of 
the tip is modeled by a square lattice with two atoms in the transverse direction. The end of 
the tip is represented by a single atom. This choice renders the system multimode, with two 
transverse modes. It is as simple to include a truly three dimensional tip, but the current will 
not be much affected by it. The schematic atomic structure of the tip is represented in Fig. [2l 
The Hamiltonian for the tip can be written as Ht = + Hq, where Hi, represents the 
bulk of the tip and Hq the tip's last atom. These two parts of the Hamiltonian are defined as 

-1 

Hi) = —V [c^(n, m)c(n — 1, m) + c(n — 1, m)c^(n, m)] (22) 

n=—oo m=l,2 
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Figure 2. (colour online) Representation of the STM tip. 



and 



-V±J^ [ct(n,l)c(n,2) + c(n,2)ct(n,l)], 



Ho = eoc^(0)c(0) -WiJ2 [cHO)c(-l, m) + c(-l, m)c\0)] , 



(23) 



(24) 



m=l,2 



where c^^ (c) are creation (annihilation) operators for fermions in the tip. Let us now consider 
the case of the bulk part of the Hamiltonian's tip, Hi,. In the case of a square lattice the wave 
function is separable and can be written as {ipi^t) = \4^i)\4^t)^ with the longitudinal part of the 
wave function given by 



-N 



i) = lim 



N^oo ^ V + 1 
n=— 1 



sin(n6'i)|n) 



and the transverse part given by 

1*^*^ = y|sin(mat)|m). 



(25) 



(26) 



m=l,2 



In Eqs. (|25l) and (|26|) . the states |n, m) = |n) |m) are position states and the numbers Qi and 
ttf are given by 



and 



1 



/ = l,2,...,iV, 



«4 = y , t = 1,2. 



(27) 



(28) 



The resolvent for the Hamihonian H}, is defined as = {E + iO~^ — Hi,) ^, where the + 
superscript denotes the retarded function. In the eigenstate basis, it has the form 



_ m,t/ 



i,t 



E — El 



(29) 



where are the eigenvalues of Hb, with Hb\i)i,t) = Ei^t\M, given by 

Ei^t = -2V cos di - 2V± cos a* , (30) 

For the calcution of the STM current we will need the surface Green's functions defined as 

Gdiaa{E) = {m,-l\Gt\-l.m), (31) 
a//,(^) = (l,-l|G+|-l,2). (32) 
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The calculation of (|3T1) and (l32l) requires the evaluation of the integral 



1^1- 

27r .,0 



2n 



which is easily done by contour integration methods 13311 . The final results are 



j = ±l 



(33) 

(34) 
(35) 



j=±i 



for /3j > 1, with /3j = {E + jV±)/{2V). In the case < 1, the Green's functions are 
obtained from ((34l) and (|35|) by removing the factor sgn(/3j), and choosing the positive sign 
for the square root of the negative argument: 



Gdiag{E) - E 
j=±l 



2V 



Goffd{E) 



j2i 

2V 



2V 
1 



1 



PI 



1 



(36) 
(37) 



^ 2V ■''2V 
i=±i 

Using Eq. (|34l) . the local density of states at the sites n = —l,m = 1,2 given as usual by 
Pb{E) = —^^GdiagiE), is depicted in Fig. [3l The multimode nature of the tip is clearly seen 
in the form of the density of states. 



1 — '"T"' — r^~\ — '~r 



1 — ' — — r^~\ — '~r 




E(eV) E(eV) 

Figure 3. Local density of states, Pb{E), at the atoms of the tip given by n = —1, m = 1,2. 
Left: y = 2 eV, VI = 2 eV. Right: F = 2 eV, VI = 1 eV. The multi-mode nature of the tip 
is clearly seen. 



3. Graphene 's local density of states 

The properties of the STM current depend on the local density of states below the tip of the 
microscope. In this section we compute the Green's function of graphene in real space, from 
which the local density of states can be obtained. In Sec. 12. 1[ the position vector r denotes 
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the position of the unit cell, which contains two atoms. The local density of states (per spin) at 
the sub-lattice A and sub-lattice B atoms of the unit cell localized in the position r is defined 
as 

p^(r, uj) = ^ImG^^(r, r, cu) , (38) 

where x = a,b, and G^xl'^, r, cj) is obtained from 

k,p 

after the usual analytical continuation ujn ^ uj + iO^ of the Matsubara Green's function. In 
the unit cell r = 0, which contains the impurity in its A site, we have simple expressions for 
the Green's functions for the A and B sites. These read 

GaaM = G'M + gM + 2hMG'{ujn) + [G'M]'TM , (40) 
= G^K) + -^[G'MTnur,) , (41) 
where G^{(jJn) is given by 

G'icUn) = -lUJnt-^ + {lUJnft-^G\uJn) , (42) 

and G^{ijJn) is defined in Eq. ([20h . As is clear from the above equations, the central quantity 
that needs to be calculated is the integrated Green's function G^{(jJn)' Since we want to 
compute the density of states for energies beyond the Dirac cone approximation we need 
to include in the density of states powers of the energy beyond the usual linear term. In a 
previous work [l35]l . we have derived an expansion for the density of states (per unit cell, per 
spin) valid for energies up to ^ 2.5 eV, reading {E = Tiw) 

Using this expression for the density of states, a close form for the retarded function 
G^{(jJn UJ + iO^) can be derived. The imaginary part of G^{uj) reads 

^G'{cj) = -^p{hw), (44) 

and the real part has the form 

m'icv) = Prihv) + P,{Iuv) In , (45) 

where Pi{x) and P2{x) are polynomial functions given by 

The energy Dc is a cut-off energy chosen as = ^/37^t'^. It is possible to derive simple 
analytical expressions for the local density of states at the unit cell r = 0, where the impurity 
is located, using the Dirac cone approximation. If the full forms, Eqs. (l44h and (1451) . of the 
Green's are used, an analytic close form is still possible, but is somewhat cumbersome. 
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For the calculation of the local density of states, the case of a vacancy and the case where 
7^ and ^ t have to be treated separately. For the vacancy, the density of states at the 
neigbouring B atom is 



Pb(0,cu) = 



1 - 



1 



+ 



1 



with 



1 + 



^3 



TT 



(48) 



(49) 



For the case of a general substituting atom, the local density of states at the impurity atom 
Pa(0, u) is obtained from the imaginary part of the Green's function which reads 

to{2t-to){t-tof 



Nc\D{ij)\^ ' ' iV,|D(a;)|2 

The imaginary part of the Green's function at the B site, next nearest neighbor to the impurity, 
reads 



1-^ 



^N,\D{uj)\^ V ujJ j ' ^^^^ 

If in Eqs. dSO]) and dSB one uses the full result for G^{uj), given by Eqs. (|44l) and (1451) . 
the resulting expressions are valid for energies up to 2.5 eV. The local density of states 
Px{r = 0, uu), X = a,b, are depicted in Fig. |4l for different choices of to and e^. 
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Figure 4. Local density of states, Px{E), for x = a, 6 at the graphene's unit cell r = 0. 
We have used t = 3 eV. Upper left: density of states of pristine graphene and at the B site 
close to a vancancy. All other panels: local density of states at A (impurity) and B (next to 
the impurity) sites for different values of the parameters to and e^. The values of to and 
corresponds to different types of impurities. 



It is well-known that the density of states due to a vacancy has a strong departure from 
the pristine value close to the Dirac point. This is due to the logarithmic singularity seen 
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in Eq. (1481) . In the case where there is an enhancement of the hopping amplitude between 
the impurity and the neighboring atoms (negative to) the local density of states retains its 
linear behavior close to the Dirac point, but its value at the A and B sublattices are different, 
as expected (Fig. |4l upper right panel). This case would mimic a boron impurity atom. 
Boron has a larger atomic radius (R 0.85 A) than carbon (R 0.7 A), and there should 
be an increase in the absolute value of the hopping amplitude when it substitutes a carbon 
in graphene. In the case of a decreasing of the electron hopping between the impurity and 
the carbon atoms (positive to) the behavior of the density of states is more interesting since 
resonances start to develop around the Dirac point (Fig. |4l two lower panels). Note that the 
density of states still goes to zero at the Dirac point. This behaviour is reminiscent of the 
vacancy, since we can picture the two resonances developed in both sides of the Dirac point 
as a splitting of the divergent peak for the vacancy due to the departure of to from its vacancy 
value to = t' This case would mimic a nitrogen atom, which has a smaller atomic radius 
(i? 0.65 A) than carbon. 

The calculation of the local density of states for finite r requires the calculation of 
the Fourier transform entering the definition in Eq. ([39l) . For Ghh{r,r,uj) we use an 
approximation for (j){k), which reads (j){k) Sao{ky — ikx)/2. Carrying out the Fourier 
transform we obtain 

Gaa{r^ r, u) = G\uj) + T{uj) [F,{uj, t)Y , (52) 
G,,{r,r,uj) = G\uj) + ^T{uj)[F,{cu,r)]\ (53) 
where r) (n = 0, 1) is defined as 



2n+l 



TTV 



~T 2 — -^ — Jn{a)H 

r"./n a"^ — vJ^ 



(54) 



with Jn{x) the Bessel function of integer order n; and kc = 2y^/ 3\/3ao), a = \(jj\r/vF, 
Ac = 3\/3ao/2, and = 3tao/2. The Cauchy principal value of the integral in Eq. ([54h is 
computed using numerical methods. 

In Fig. [5l we plot the local density of states at both sub-lattices A and B. The typical 
oscillations due to the presence of the impurity are present. Note that close to the impurity 
the B sublattice density of states presents higher harmonics as function of r; these are due to 
the cut-off momentum kc. On the other hand, the large wavelength oscillations are the 2qF 
Friedel oscillations: from Fig. [5]the wavelength is about A 28ao; on the other hand, for the 
energy huj = 0.5 eV the Fermi momentum is qp = l/(9ao), implying A t^/qf — 28ao. 
At large values of r the two density of states are out of phase by a factor of tt. Therefore, 
when we average over the unit cell the result is essentially the pristine density of states. This 
result has strong consequences for STM experiments. If the STM experiment lacks atomic 
resolution at the A and B sublattices level, the experimental data will show a very faint trace 
of the 2qF Friedel oscillations. This seems to be the case in the experiments reported in Ref. 
|^8J. Closer to the impurity there is a strong departure from the asymptotic behavior [|9j] 

p(r) oc ^8m{2ru/vF) , (55) 
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Figure 5. Local density of states in sublattice A and B as function of r > 0. The dashed- 
dotted line is the density of states of pristine graphene at the energy huj = 0.5 eV. The point 
r = is excluded. The small wavelength oscillations are due to the cut-off momentum kc and 
the large ones are the 2qF Friedel oscillations. 



and the A and B LDOS behave quite differently. 

STM measurements of a material surface are ideal for studying real-space local features 
with atomic resolution. In particular, real space modulations of the STM intensity can be 
observed when impurities are present at the surface of a given material. In general, the 
impurities lead to elastic scattering between the momentum and —Qf, which is the most 
efficient process due to phase space restrictions, leading to 2qF Friedel oscillations. In the 
case of graphene, the chiral nature of its electronic spectrum changes this general behavior. 
The Fermi surface has disconnected pieces at different points of the Brillouin zone - the K 
and K' points. The Fermi surface consists of circumferences of radius Qf around each K 
and K' points. The scattering process is then characterized by two channels: an intra-cone 
scattering (within the same K or points) of momentum change 2qF, and an inter-cone 
scattering (between the K and K') points. 

A Fourier transform of the real-space STM-intensity currents, proportional to the LDOS, 
will produce bright spots at the momentum values seen in the real space modulations of 
the LDOS. When impurities are present, the momentum values characterizing the real space 
modulation are related to the momentum change associated with a given scattering process. 
In Ref. [28J it was found that the intra-cone scattering, which would give rise to a bright spot 
of radius 2qF, was absent in the momentum map of the density of states obtained by a Fourier 
transform of their STM data. 

In what follows, we give Fourier transforms of the local density of states of graphene for 
the different types of impurities discussed previously in the text. We note that our derivation 
of the Fourier transform of LDOS uses the full Green's functions for the calculation of the 
LDOS, and therefore no approximation has been made in the calculation. The full real-space 
map of density of states can be obtained numerically from Eqs. ([38]) and ([39l) . using the exact 
expressions Eqs. ([T4h-([2T1). as was done in Ref. ifTOll . We now perform a Fourier transform of 
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Figure 6. Fourier transform of the local density of states. The first row is the density of states 
in sublattice A (pa), the second row is the density of states in sublattice B (p^), and the third 
row is the sum of the two, pa + Pb- Column (a) is for the case of a vacancy, and co = 0.3 
eV. Column (b): to = -1 eV, = 1 eV, co = 0.3 eV. Column (c): to = 1 eV, = -1 eV, 
= -1.5 eV. Column (d): to = 2 eV, = -1 eV, u; = -0.3 eV. 

the density of states, 

u;)=Yl e-^^-V^(r, u) , (56) 

r 

for each sublattice x = a,b. The results are shown in Fig. [6l where the three rows correspond 
to pa{k,(jj), pb{k,(jj), and the sum pa{k,(jj) + pb{k,(jj), from top to bottom; and the four 
columns correspond to four types of impurities discussed in Fig. |4l We recall that positive to 
reduces the hopping from the impurity site to its nearest neighbor carbon atoms - the particular 
case of to = t represents a vacancy - and negative to increases the hopping of the electrons 
from the carbon atoms to the impurity site. This latter case would correspond to an atom with a 
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Figure 7. Fourier transform of the local density of states for to = 2 eV, = — 1 eV. The 
first row is the density of states in sublattice A(pa), the second row is the density of states in 
sublattice B (p^), and the third row is the sum of the two, pa + Pb- Column (a) is for co = 0.05 
eV. Column (b): u; = 0.15 eV. Column (c): u; = 0.3 eV. Column (d): u; = 0.5 eV. 



radius larger than carbon, such as boron, leading to an increase of the hopping relatively to the 
hopping t between nearest neighbor carbons. The column (a) of Fig. [prefers to a vacancy. In 
this case, it is clear that a 2qF circumference is seen around the k = (0, 0) point for the pa and 
P5 plots, consistent with the modulations shown in Fig. O However, the intensity at the 2qF 
circle around k = (0, 0) is suppressed when looking at the pa + Pb plot. Features at six spots at 
a distance \K\ around k = (0, 0) are also seen, these correspond to the K and points at the 
corners of the first Brillouin zone, and represent inter-cone scattering. Additionally, six bright 
spots at distance are also present. These vectors with modulus are reciprocal 

lattice vectors G. For a pristine material the local density of states has the periodicity of the 
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underlying lattice, that is, p(r) = p{r + R), and therefore, the Fourier transform of p(r) must 
show the same intensity at fe = (0, 0) and k = G. While the presence of an impurity breaks 
the periodicity of the real-space lattice, the reason for k = (0, 0) and k = G being different 
is the fact that the impurity considered here is not on a whole unit cell, but on only one of the 
sites of the unit cell. If there were only one site per unit cell and a single short-range impurity, 
then the periodicity in k-space would be maintained. In the extreme case, if there was a line 
of impurities of the vacancy type, this would correspond to cutting the system into half, and 
k = (0, 0) and k = G would still be the same. Note that in all our cases pa is still periodic 
in /c-space. The impurity on a site A that we consider here, introduces structure within the 
R = unit cell, therefore, in k-space there is information going beyond the first Brillouin 
zone. Mathematically this appears because the B-site is located at R + Ss, where R denotes 
the position of the A sites, and 63 is not a lattice vector of triangular Bravais lattice (see Fig. 
[B. When we take the Fourier transform, the periodicity in fc-space does not happen at fc = G 
anymore, but at at larger k. All figures show this signature, which is rather clear in column 
(c) of Fig. [6l since the Fermi surface energy has been chosen as large as = — 1.5 eV. 

Inter-cone scatterings, represented by the region around the K and points, are highly 
angular-dependent [|29l[3^, as can be seen particularly in the plots. While the scattering 
around the G vectors do not have such a strong angular dependence, some trigonal warping 
is observed in p^. The scatterings around k = (0,0) are rotationally symmetric. Fig. |7] 
corresponds to the case to = 2 eV and = — 1 eV (case (c) of Fig. [6]), showing how the 
/c-space LDOS map evolves with increasing energy. 



4. STM current 



In this Section we present calculations of the tunneling current between the STM tip and 
graphene, when the tip is close to an impurity atom. We model the tip by the multimode 
tight-binding model, as described in Sec. 12. 2[ This choice departs from the more simplified 
approach where the tip is modeled by a one dimensional system f36U37L 

There is a number of ways one can use to describe the tunneling of the electrons between 
the STM tip and graphene. Here we assume that the coupling is made directly either to the 
impurity atom or to the next neighbor carbon atom. This choice corresponds to probing the 
local electronic properties at or around the impurity. More general types of coupling are easily 
included in the formalism. We write this coupling as 

Ht = -W2[c\0)d{0) + rft(0)c(0)] , (57) 

where the operator d{0) can represent either an electron at the impurity atom in the A sub- 
lattice or at the carbon atom in the B sub-lattice. 

Since the Hamiltonian of the problem is bilinear we can write it in matrix form as 



H, Vr 



H 



R 



(58) 



Vr 
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where the matrices Vl and Vr represent the couphng of the last atom in the tip of the STM 
microscope to the bulk of the tip and to graphene, respectively, and and Hg stand for the 
bulk Hamiltonians of the tip and of graphene, respectively. Hg also includes the impurity 
terms Eqs. ^ and (|3]). The matrix H is of infinite dimension due to and Hg. The matrices 
vl and have the explicit form 

vl = [0, -H^i, -W^] , vl = [-W2^ 0] , (59) 

where represents an infinite dimensional null row vector. 

The tunneling is a local property, controlled by the coupling of the last atom of the tip 
to the bulk atoms and to graphene. Since we want to compute local quantities, this is best 
accomplished using Green's functions in real space. The full Green's function of the system 
is defined by 

{lE + iO^ - H)G^ = 1, (60) 
where 1 is the identity matrix. The matrix form of the Green's function is 





Gbb 


Gbo 


Gbg 






Gob 


Goo 


Gog 


(61) 




Ggb 


GgO 


^99 





The quantity of interest is Goo, which can be shown to have the form 

G^, = {E + zO+ - 60 - S+ - S+) , (62) 
where the matrices and are the self energies and have the form 

E+ = 2W^iGa^a, + Goffd) , S+ = W^G^, , (63) 

where G^^ is the surface Green's function of the Hamiltonian Hg at the impurity unit cell 
(x = a, 6), respectively. Note that the quantity G^^ is computed using Eq. (1391) and setting 
r = 0. 

The study of non-equilibrium transport is done using the non-equilibrium Green's 
function method, or Keldysh formalism. This method is particularly suited to study the regime 
where the system has a strong departure from equilibrium, such as when the bias potential on 
the STM tip, V^, is large. In this work, we consider, however, that the system is in the steady 
state. Since the seminal paper of Caroli et ah on non-equilibrium quantum transport [38J, that 
the method of non-equilibrium Green's functions started to be generalized to the calculation of 
transport quantities of nanostructures. There are many places where one can find a description 
of the method OH HOI, but a recent and elegant one was introduced in the context of transport 
through systems having bound states, showing that the problem can be reduced to the solution 
of an equation similar to a quantum Langevin equationp4P|. The general idea of this method 
is that two perfect leads are coupled to the system, which is usually called the device. In our 
case the device is defined by the last atom of the tip of the microscope. The Green's function 
of the device has to be computed in the presence of the bulk of the tip and of graphene. This 
corresponds to our Gjo Green's function. Besides the Green's function we need the effective 
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coupling between the last atom of the tip and the bulk atoms as well as that to the graphene 
atoms, which are determined in terms of the self-energies 

^L/R=^{^t/R-^-L/R)- (64) 

Therefore the effective coupling T^/ji depends on the surface Green's function of the tip and of 
graphene. According to the general theory, the two systems (bulk of the tip and graphene) are 
in thermal equilibrium at temperatures Tl/r and chemical potential fiL/R and are connected 
to the system at some time to- The total current through the device is then given by 

2e f"^ 

J=-^ J dET{E) [f{E, fiL. Tl) - f{E, fiR, Tr)] , (65) 

where the factor of 2 is due to the spin degrees of freedom, f{x) is the Fermi-Dirac distribution 
and the transmission T{E) is given by 

T{E) = i7r'TL\G^,\'TR. (66) 

In Fig. [8] we depict T{E) in different cases. Note the asymmetry of the density of states 



' I ' I ' I ' ; ' I ' 




E (eV) E (eV) 

Figure 8. Transmission probability T{E). The parameters used are (all in electron-volt): 
V = 2,V± = h 1^1=0.9, W2 = 0.2, and eo = 0.2. The resonances seen in the LDOS in Fig. 
Hlshow up in the transmission function. The values of to and are the same used in Fig. |4l 
and the four panels here correspond to the same ones in that figure. 



which is exhibited even by the pristine case (Fig. [8l upper left panel). This asymmetry has 
a two fold nature: (i) it comes from the fact the bulk of the tip has two transverse atoms but 
the tip has only one; (ii) the fact that the atom at the tip has a different local energy from 
those in the bulk. This asymmetry carries on to the disordered cases. Additionaly, for the 
disordered cases the resonances seen in the local density of states has a strong impact on the 
transition probability T{E), leading to open transport channels with large values of T{E). 
This is specially true for the vacancy and for the weakly coupled impurity case, that is, when 
hopping between impurity and carbon atoms is suppressed in relation to the hopping between 
carbon atoms, as expected for nitrogen substitution. 

Since we want to probe the properties of the STM current at zero doping we choose 
jJiL = eV/2 and jir = —eV/2. Also Tl = Tr = 0. This renders the calculation of the current 
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V^(volt) V^(volt) 



Figure 9. STM current J. The parameters used are (all in electron-volt): V = 2, V± = 1, 
Wi=0.9, W2 = 0.2, and eo = 0.2. The values of to and are the same used in Fig. IH The 
transmission function was computed at finite bias; the zero bias case is given in Fig. [S] The 
four panels here correspond to the same ones in that figure [H 



to a simple one-dimensional integral of T{E) over the energy. The form of the current will 
reflect the properties of T{E) as function of energy, and, as we have seen, these are markedly 
different for the different cases, depending strongly on the value and sign of to. The presence 
of resonances in T{E) leads to steps in the STM current. This is seen in Fig. [9] for the cases 
of the vacancy and to the case of weak coupling (positive to) between the impurity and the 
neighboring carbon atoms. 

To fully characterize the STM current, another important quantity is the shot noise P2I]. 
For interacting systems, this quantity contains information on the nature of the quasi-particles, 
including for example the possible existence of quasi-particles with fractional charge. In 
disordered systems with no interactions, information on transport through open channels can 
be obtained. For non-interacting electrons at zero temperature the shot noise is defined as HSll 

S=— dE T{E) [1 - T{E)] . (67) 

The relevant quantity is not S directly but the Fano factor ll42l defined as F = ^. When the 
transmission T{E) is strongly reduced we have F ^ 1, and the noise is said to be Poissonian. 
On the other hand, if the system has a finite density of open channels, T{E) 1, we have 
F < 1 due to [1 — T{E)] <C 1. The resonances that we obtain in T{E) play a role in the 
resulting form of F. They lead to an enhancement of the current, and to a significant decrease 
of the Fano factor due to the opening of a transport channel. Note that the opening of the 
channels is a consequence of the local disorder induced by the impurity. 

5. Conclusions 

In this paper we have studied the STM currents through locally disordered graphene. We have 
considered a tip with transverse modes. Although the tip is strictly quasi-one-dimensional it 
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V, (volt) (volt) 



Figure 10. Fano factor F. The parameters used are (all in electron-volt): V = 2, V± = 1, 
Wi=0.9, W2 = 0.2, and eo = 0.2. The values of to and are the same used in Fig. IH The 
four panels here correspond to the four panels of the current given in Fig. [51 

Still departs from the widely used model of a strictly one-dimensional model. Generalizing 
now the calculations to a truly three-dimensional tip is reasonably straightforward. The 
modifications would require introducing a three dimensional square lattice for describing the 
bulk of the tip, and a decreasing number of atoms for each transverse plane to describe the 
part of the tip in contact with graphene. This last part would lead to the most significant 
change in the calculation, since the device would not be a single atom as in our calculations 
here, but would be represented by a finite number of them, and therefore the Green's function 
for the device would be a matrix instead of a c-number. Nevertheless, as long as we take the 
dispersion of the electrons in the tip to have large bandwidth, the current should not depend 
much on the local density of states of the tip, since this would be essentially constant. This 
corresponds to the usual wide band limit. 

We have also seen that tunneling through either impurity atoms, or their neighboring 
carbon atoms, depends on the local density of states of graphene. For certain circunstances 
- vacancy or weakly coupled impurities - there is a development of resonances at or close to 
the Dirac point. These resonances lead to a strong enhancement of the tunnelling probability 
which appear as steps in the tunneling current. It is conceivable that graphene could be locally 
modified in order to take advantage of these strong resonances developed close to the Dirac 
point. Clearly the substituting atoms would also locally distort the graphene lattice, an effect 
not included in our description. How much the STM current would depart from the values 
computed here would depend on the change in the values of the hopping parameter, and on 
additional features on the density of states due to the disorder. If future research will pursue 
the route of modifying graphene locally, our results will be important for the characterization 
of the surface. Even in the present state of affairs, our results could be used to interpret STM 
current due to local impurities. 
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